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ABSTRACT 



For pulsar projects it is often necessary to predict the pulse phase in advance, 
for example, when preparing for new observations. Interpolation of the pulse phase 
between existing measurements is also often required, for example, when folding X- 
ray or gamma-ray observations according to the radio pulse phase. Until now these 
procedures have been done using various ad hoc methods. The purpose of this paper 
is to show how to interpolate or predict the pulse phase optimally using statistical 
models of the various noise processes and the phase measurement uncertainty. 

Key words: pulsars: general - methods: data analysis 



1 INTRODUCTION 

Pulsars are rapidly rotating neutron stars. Each rotation of 
a radio pulsar leads to the characteristic pulse of radiation 
that can be detected using a radio telescope. Because the ro- 
tation of pulsars is exceptionally stable, the time of arrival 
(ToA) of each pulse can be precisely predicted from a timing 
model that includes the motions of the pulsar and the Earth 
and the various propagation effects between the pulsar and 
the Earth. Comparing the actual pulse arrival time with the 
prediction of a timing model, whose parameters are adjusted 
using a least-squares procedure, is the essence of pulsar tim- 
ing. The timing model and least squares proce dure are im- 
plemented in software pack ages, e.g., tempo2 (|Hobbs et al.l 
2006; Ed wards et alj|200d 1. A preliminary timing model is 
obtained with the discovery of the pulsar. When a series 
of observations have been made over a period of years, the 
timing model can be refined with great precision. The devia- 
tions between the actual and the predicted arrival times are 
known as "timing residuals." Ideally, these residuals would 
simply reflect the ToA measurement uncertainty, but in fact 
they will include various un-modeled effects. Residuals often 
exhibit slow variations which appear to be a red stochastic 
process, which is referred to as "timing noise" . Timing noise 
may be caused by many phenomena including those intrinsic 
to the pulsar (e.g., Hobbs et al. 2010, Lyne et al. 2010 and 
references therein), effects in the interstellar medium (e.g., 
You et al. 2008), the interplanetary medium (You et al. 2012, 
in press) or even irregularities in terrestrial time standards 
(Guinot & Petit 1991, Petit & Tavella 1996, Rodin 2008 and 
Hobbs et al. 2010). 



Being able to remove the timing noise is impor- 
tant for numerous applications. For instance, in order for 
iHobbs et aD (|2004h to measure accurately pulsar proper mo- 
tions, it was necessary to remove timing noise without affect- 
ing any signal with a periodicity close to 1 yr. This was done 
using the fitwaves algorithm (described in the Appendix 
of Hobbs et al. 2004) which fits a sequence of harmonically 
related sinusoids to the data, where the period of the high- 
est fre quency wave was constrained to be > 1.5 yr. Recent 
work (|Coles et alj|201ll ) showed that improved parameter 
measurements could be made by obtaining and using a sim- 
ple, analytic model of the power spectrum of the timing 
noise in order to calculate the covariance matrices of the 
red and white noise processes. With the covariance matri- 
ces one can find a linear transformation that whitens the 
observations. If this transformation is applied to both the 
observations and the timing model, the least squares prob- 
lem is greatly simplified. With the covariance matrices the 
red timing noise can also be interpolated or predicted us- 
ing a maximum likelihood estimator. It is this interpolation 
and prediction scheme that we will discuss in the following 
sections. 

In order to obtain the characteristic pulse profile for a 
pulsar, it is necessary to sum many individual pulses dur- 
ing the observation. For radio observations this is often car- 
ried out "online" by folding the incoming data stream at 
the known topocentric period of the pulsar. Over short data 
spans, this can be done with sufficient accuracy by assum- 
ing that the timing model is perfect. However, if a pulse 
profile is obtained from a data set that spans many months 
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or years then it is often necessary to model the timing noise. 
For instance, obtaining a pulse profile using the Large Area 
Telescope (LAT) on the Fermi Gamma-ray Space Telescope 
requires the gamma-ray photons to be added together ac- 
cording a timing m odel that must be valid over many years 
l|Abdo et al.ll2010f ). This is commonly done using a timing 
model from observations of the pulsar obtained with a radio 
telescope and then "folding" the gamma-ray photons into 
the predicted phase bin according to this model. The radio 
model gives timing resid uals that may be "white ned" using 
the fitwaves procedure l|Weltevrede et alj |2010l. However, 
as we demonstrate in subsequent sections, the fitwaves pro- 
cedure is not theoretically optimal in the sense that a max- 
imum likelihood estimator is optimal, nor can it be extrap- 
olated past the end of the radio data (or back-extrapolated 
before the start of the radio data). Fitwaves also requires 
an arbitrary choice of the number of harmonics to include 
in the fitting procedure. Too few harmonics implies that not 
all the features in the timing residuals are completely mod- 
elled. Too many harmonics leads to an unphysical model for 
the timing residuals (particularly when large gaps exist in 
the data). 

Prediction of the pulse phase into the future is required 
for real-time folding of pulsar data (a very widely used tech- 
nique). This is typically carried out by extrapolating the 
pulsar timing model into the future and ignoring the ex- 
pected variations due to the pulsar timing noise (the extrap- 
olated timing model is typically described using a Chebyshev 
polynomial ex pansion for efficien t implementation in digital 
hardware; see lHobbs et al. I l2006l for details). 

Another example of a practical use of the requirement 
to predict pulse phase into the future is in order to help 
navigate space probes in the Solar System (ISheikhl 120051 ; 
ISheikh fc Pmesll2006l : iBernhardt et al.ll20ld . l201lf ). Compar- 
ison of pulse arrival times at the satellite and the expected 
pulse arrival times given a pulsar timing model can be used 
to determine the position of the satellite. However, it is ex- 
pected that any such system will only have the sensitivity 
to observe pulsars with high X-ray flux density; such pulsars 
are typically young and have residuals that exhibit signifi- 
cant timing noise (e.g., Gavriil & Kaspi 2002). Optimal ex- 
trapolation of the timing residuals is therefore essential for 
such navigation. 

In this paper, we describe how to use the covariance 
matrices of the red and white noise processes to both in- 
terpolate and extrapolate pulsar timing residuals, without 
the need for fitwaves. In Section 2, we outline the theory. 
In Section 3, we carry out simulations to demonstrate the 
effectiveness of our technique. In section 4, we compare the 
new technique with fitwaves in timing noise modelling of 
radio data from observations of the Vela pulsar, B0833— 45 
for folding of Fermi gamma ray observations and outline 
the advantage of the new algorithm in space-craft naviga- 
tion. In Section 5, we apply our algorithm to a simulated 
data set that includes a glitch event, use simulations to es- 
tablish the sensitivity of the method to errors in param- 
eterising the noise and use real data from observations of 
PSR J1939+2134 to compare the new algorithm with and 
without quadratic removal. We conclude the paper in Sec- 
tion 6. 



2 THE MAXIMUM LIKELIHOOD FILTER 

The theory of optimal interpolation and extrapolat ion o f 
stationary time series was first discussed by Wiener (|l964l ) 
and the corresponding algorithm has become known as a 
Wiener filter. However normal Wiener filters are not appli- 
cable to pulsar timing residuals because the measurement 
uncertainty is not stationary and the sampling times are 
irregular. However, Wiener filters are maximum likelihood 
estimators (MLE) and we can obtain an MLE for pulsar tim- 
ing residuals directly as shown below if we know the covari- 
ance of the red and white noise processes. This is effectively 
a generalised Wiener filter. 

Here we wish to estimate a red noise process (s), from 
irregularly spaced observations (o) of the red noise plus mea- 
surement error (n), i.e. o = s + n, when both processes are 
gaussian. The likelihood function, i.e. the logarithm of the 
probability density, for N observations (neglecting the nor- 
malizing constants) is given by 

s T C^~ 1 s + n T Cn 1 n = s T C^ 1 s + (o — s) T C^ 1 (o — s) (1) 

Here C n is the covariance matrix of the measurement error 
and C s is the covariance matrix of the red noise. To find 
the MLE estimator for s given observations o, we set the 
gradient of the likelihood function with respect to s to zero 
and solve the resulting linear system, 

2C~ 1 s-2C~ 1 (o-s) = (2) 

(C" 1 + C-> = C-V (3) 

C n is a diagonal matrix so it is trivial to invert, but the 
inversion of C s can be unstable so the Cholesky decomposi- 
tion should be used if it is necessary to invert it. Equation 

3 can be solved directly or it can be simplified first by pre- 
multiplying both sides of C n or C s . We pre- multiply both 
sides by C s to avoid the need for inverting C s . 

This solution provides an estimate of the timing noise 
s at every observed sample. However the MLE is not as 
restrictive. We can solve for the signal that maximizes the 
likelihood function, even at locations where we do not have 
a sample o. In this case the likelihood function will include 
s terms that are not matched by o terms. Thus in equation 
2 the term C„ 1 (o — s) will include only the subset of s 
which matches o. We simplified the system by augmenting 
the o vector with terms of value zero terms matching the 
unmatched s terms. We also augmented C^ 1 with terms 
of value zero along the diagonal, corresponding to the aug- 
mented terms of o. With these extensions to o and C„ 1 the 
system can be reduced to the form in Equation 3 and solved 
as discussed earlier. 



3 TESTING THE ALGORITHM 

In order to test our algorithm we first simulated data sets 
that contain both white and red noise. In order to obtain 
realistic sampling and ToA uncertainties we used observa- 
tions of PSR J1713+0747 obtained using th e Parkes radio 
telesc ope from MJD 49421.9 to MJD 54546.8 l|Verbiest et all 
2009) as an example of a "typical millisecond pulsar" data 
set. With this sampling, we simulated a data set of timing 
residuals with red noise and added white noise to these gen- 
erated observations. We carried out simulations with three 
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Figure 1. Power spectral densities for the three types of "timing 
noise". The solid line is defined by A = 7.6 X f0 _30 ,a = 4.3333, 
f c = O.lyr -1 and f a = lyr -1 . The dashed line with A = 7.6 X 
fCT 30 , a = 4.3333, / c = O.lSyr" 1 and f a = fyr" 1 . The dotted 
line with A = 7.6 X 10- 30 ,a = 4.3333, f c = O^Syr" 1 and f a = 
fyr - 1 . The vertical dot-dashed line is at a frequency of f/T span . 




5000 10* 1.5xl0 4 

Time (days) 

Figure 3. 50 realisations of each type of timing noise shown in 
Figure f. The top panel contains the results for the strongest 
(defined as the highest amplitude noise in Figure f) red noise, 
the middle panel contains the intermediate-level noise and the 
lowest panel contains the weakest red noise. The rms of these re- 
alizations for each sample is shown as the dashed lines. The rms 
difference between the known simulated noise and the interpo- 
lated/extrapolated predictions are shown as the solid line. 




Time (days) 

Figure 2. The upper panel shows one realisation of simulated 
timing residuals (with same error bars as real data). The irregu- 
lar solid line shows the extrapolation of the data for 5000 d before 
and after the end of the actual sampling. The dotted line indi- 
cates the prediction of the timing noise using our algorithm. The 
lower panel shows the difference between the simulation and the 
prediction with the error bars overlaid where sampled. 



different forms of timing noise. In all cases the white noise 
level remained the same. The simulated spectrum model has 
the form P(f) = A/[(f c /f ) 2 + {f / f ) 2 ] {a/2) . A is the am- 
plitude of the red timing noise, f c is a corner frequency that 
ensures that the red noise model turns over at low frequen- 
cies, f = lyr -1 and a is the spectral exponent of the red 
noise. The spectral densities of the three timing noise types 
are shown in Figure [1] We define these three noise types 
as having "strongest", "intermediate-level" and "weakest" 
red noise depending on the amplitude of the red noise pro- 
cess. Note we choose the spectral exponent of a = 4.3333 
to represent the red noise as this is expected if the timing 
residuals were dominated by the effect of a stochastic back- 
ground of gravitational waves (e.g., Sesana, Vecchio & Cola- 



cmo Similarly steep spectral exponents are also typ- 

ical of other forms of timing noise such as residuals caused 
by unexplained rotational irregularities , glitch recovery or 
unmo delled interstellar medium effects (jCordes fc Shannon! 
l20ld ). The cutoff frequency was chosen to lie close to 1/Tspan 
(where T span is the data span of the actual observations), 
shown as the vertical dot-dashed line. 

For each of the three spectral models we computed a 
single daily-sampled realization of the timing noise. In or- 
der to test our extrapolation procedure, each simulation had 
a data span of the actual observations, but with an extra 
5,000 days added to the start and to the end of the data 
span. We require multiple realizations of this red noise. To 
do this we created one data set that was 50 times longer 
than the required data length. We created one long data set 
so that we could simulate the low frequencies (/ < 1/Ttot 
where T to t is the required data span) with a simple discrete 
Fourier transform. This data set can then be subdivided 
into blocks of length T to t to give different realisations of the 
random process. For the central region that corresponds to 
the actual observations, we interpolated the timing noise to 
the actual sample times using a constrained cubic interpo- 
lator. We subsequently added white noise to these interpo- 
lated samples at the level given by the ToA uncertainties for 
each observation. We also regularly sample the data every 
Tapan/100 ~ 50 d over the entire data span, but do not add 
white noise to these data as the regularly spaced samples 
are used to provide a comparison between the actual timing 
noise and the interpolated timing noise. In order to simulate 
the effect of fitting for the pulsar timing model, we fitted, 
and subsequently removed, a quadratic polynomial to the 
data that represented the actual observations. 

We use our algorithm to determine the signal using 

1 Note that the characteristic strain of a gravitational wave back- 
ground has the form h c (f) = A g f a s where a g = —2/3. The 
power spectral density of the induced timing residuals has the 
form P = Af 2a s- 3 . 
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the same gridding as the regularly sampled data described 
above. In this way we could compare the results of our in- 
terpolation, extrapolation and back-extrapolation with the 
known red noise that was simulated. For these simulations, 
we did not obtain the covariance functions from the data; 
instead we used the known analytical power spectra to com- 
pute the covariance functions and the corresponding matri- 
ces. The results for the intermediate timing noise signal (the 
middle spectrum in Figure [TJ are shown in Figure [5] The top 
panel shows the simulated residuals for one realization of the 
red noise (shown with error bars). The timing noise realiza- 
tion is shown as the irregular solid line. Note that the timing 
noise diverges from the timing model after the range of the 
quadratic fit. Our algorithm leads to the interpolation and 
extrapolation that is shown as a dotted line. The difference 
between the simulated red noise and the predicted function 
is shown in the lower panel of Figure [2] Error bars repre- 
senting the white noise level are overlaid, at the locations 
of the actual sampled data points. During the range of the 
observations, the interpolation works well. Discrepancies be- 
tween the interpolation function and the simulated red noise 
level are smaller than the white noise level. Outside of the 
observing span, the extrapolated function converges to the 
prediction of the timing model (i.e., predicted residuals of 
zero). The discrepancies between the extrapolated function 
and the simulated data grow as the extrapolation time in- 
creases. The rms of the discrepancies between the model 
and the known signal measured over all 50 realizations of 
all three types of timing noise is shown as the solid line in 
Figure [3] In all cases the rms of the discrepancies converge 
to the rms of the timing noise (shown as the dashed line in 
Figure O). 



4 APPLICATIONS 

4.1 Timing noise modelling and removal 

In order to test our algorithm using a real pulsar data set 
that exhibited a significant amount of timing noise we made 
use of radio observations of the Vela pulsar (PSR B083 3— 45) 
from the Parkes 64-m radio telescope (|Weltevrede et all 
120101 ). These data were taken between MJD 54689 and 55758 
mainly to support the Fermi space-craft mission. In order to 
determine the covariance function o f the red noi s e, we used 
the Cholesky algorithm described in IColes et all l|201ll ) and 
implemented as the spectralModel plugin to Tempo2. 
We used this resulting covariance function to obtain pre- 
cise astrometric and pulse parameters for the Vela pulsar 
and formed the post-fit timing residuals. These are shown 
in the top panels of Figure|4^ and[4p (the ToA uncertainties 
are smaller than the size of the symbol). In common with 
most analyses of these data, we used the fitwaves algorithm 
(with 10 harmonics) to model the timing noise. This model 
is shown as the dotted line in the top panel of Figure [4^, 
and these post-fit residuals, once this model is removed, are 
shown in the bottom panel of Figure [4^,. Note that the data 
are modelled well during the observations, but the fitwaves 
function diverges past the end of the timing residuals and 
therefore cannot be used to extrapolate the behaviour. 

For comparison, we have applied our procedure and pro- 
vided a model for the timing noise with a grid spacing of 



13 d. The model is subsequently interpolated, using a lin- 
ear interpolation, to determine the pulse phase at the exact 
time of each observation. Our model is shown as the dotted 
line in the top panel of Figure |4p. We note that this method 
models the data well during the observations and converges 
to the predictions of the timing model about one hundred 
days after the end of the data span. The post-fit residuals 
clearly have a lower rms timing residual than the post-fit 
residuals obtained with the fitwaves algorithm. 

A significant problem with the fitwaves algorithm is 
shown in the top panel of Figure[4ji. For this Figure we have 
removed a few hundred days of radio observations near the 
start of the data span (these "removed" data points are indi- 
cated in the Figure using small triangle symbols, but are not 
included in the analysis). The fitwaves model (shown as the 
dotted line) interpolates poorly during the data gap. This 
can lead to significant problems with folding, e.g., gamma- 
ray photons that were observed during this interval. In con- 
trast, our algorithm (Figure 3J1) provides an interpolation 
that more closely represents the original observations dur- 
ing the data gap. 



4.2 Spacecraft navigation 

It may be possible to use observations of pulsars to assist 
with navigation of space probes throughout the Solar Sys- 
tem. This method relies on the ability to predict the arrival 
times of pulses precisely and accurately. However, previous 
work JSheikhl [20051 ; ISheikh fc Pines! 120061 ; iBernhardt et al.l 
|2010| . 1201 ll ) has assumed that the pulse arrival times can 
be modelled using a standard pulsar timing model. Even 
though such a timing model can be "whitened" using ex- 
isting techniques (i.e., the fitwaves procedure), it is likely 
that, for any navigational purpose, it will be necessary to 
extrapolate the timing model. As shown above, the fit- 
waves procedure cannot be used for such extrapolation. In 
Figures [4^ and [4j we have shortened the actual radio ob- 
servations of the Vela pulsar by 220 days and re-fitted for 
the pulsar's pulse frequency and its first time derivative. We 
then use fitwaves (Figure [4^) and our new algorithm (Fig- 
ure |4f) to predict the timing residuals and compare with 
the actual data. The actual data included in the analysis 
are shown as circles. The data that have been removed are 
indicated using triangle symbols. Clearly the extrapolation 
is not perfect with our method. After only a few tens of 
days the prediction has diverged from the actual data by an 
amount larger than the white noise error bars. However, in 
contrast to the fitwaves procedure, the prediction is never 
worse than the rms of the actual timing residuals. 

The error in determining the space-probe position is 
related to the error in predicting the pulse arrival times. 
However, without prior information, positional determina- 
tion can only be obtained by comparing the actual and pre- 
dicted arrival times for three or more pulsars. If the pre- 
dictions for one, or more, of these pulsars is incorrect (for 
instance, due to poor extrapolation) then it may be impos- 
sible to obtain a solution for the position of the space-probe. 
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Figure 4. Models of the timing residuals for the Vela pulsar. The left-hand panels (a,c,e) show models of the residuals using the 
FITWAVES procedure. The top part of each panel contains the original timing residuals and the model (dotted line). The lower part of 
each panel shows the difference between the actual residuals and the model. The right-hand panels (b,d,f) show the same, but using the 
interpolation/extrapolation procedure described in this paper. The panels (c,d) have data removed before the fitting, interpolation and 
extrapolation. These "removed" observations are shown with the small triangle symbols. The panels (e,f ) have the most recent 220 d of 
data have removed from the analysis. 
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Day 

Figure 5. Interpolation of timing residuals that include a glitch 
event. The upper panel shows the timing residuals with the in- 
terpolated function shown as a solid line. The lower panel shows 
the difference between the actual data and the interpolation. The 
error bars are defined by the 100 ns of white noise simulated for 
this data set. 
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Error in model parameter 

Figure 6. Sensitivity of the measurements to changes in the noise 
model. The dotted, dashed and solid lines show the sensitivity to 
changes in the parameterisation of the spectral component, a, 
corner frequency, f c and the white noise, a^, respectively. These 
are represented on the horizontal axis by a — ao, log 1 Q(/ c // c o) 
and log 1 o(o"2j/o'^ j0 ) where aoi fcO and ct^q are the known values 
for these parameters. 



5 DISCUSSION 
5.1 Glitch events 

The algorithm presented here is correct if the covariance 
matrix exists. It does not have to be stationary (or wide- 
sense stationary). However the scheme for estim ating the 
covariance matrix presented bv lColes et al.l lj2011T j assumes 
that the red noise is stationary. If the noise is non-stationary, 
for example if small glitch events exist in the residuals, then 
this information should be included in the timing model to 
the extent possible. However, in some cases, the glitch event 
may not be easily identifiable. Here we consider whether 
the covariance matrix can deal with weakly non-stationary 
observations. 

We have simulated a data set with regularly spaced ob- 
servations at 14 d intervals and white noise with an rms tim- 
ing residuals of 100 ns. The exact amount of white noise is 
not important; this value was chosen so that the effect of the 
glitch is large compared with the white noise level. We have 
added a moderate red noise signal to simulate the effect of 
timing noise and included a small glitch (corresponding to 
a frequency change of 1 x 10~ 12 Hz) exactly at the centre of 
the data span. We have formed timing residuals, obtained a 
spectral model of the noise and used our algorithm to inter- 
polate the data points. The result is shown in Figure [5] The 
difference between the interpolated timing residuals and the 
simulated data is less than the white noise level, demon- 
strating that the interpolation procedure correctly models 
the timing residuals. 



5.2 Sensitivity to errors in the parameterisation 
of the noise 

Our algorithm allows pulsar timing residuals that exhibit 
red noise to be modelled by interpolation or extrapolation. 
In contrast to previous methods, such as using the fitwaves 
algorithm or fitting high-order polynomial expansions, the 



extrapolated model converges on the prediction of the pulsar 
timing model. 

Our MLE method is optimal only if the statistical prop- 
erties of the noise is known, so it is important to establish 
its sensitivity to errors in the parameterisation. We have re- 
peated the simulations described in Section 3, but here we 
used incorrect estimates of the spectral model and the white 
noise covariance. The timing noise simulated corresponds to 
the lowest spectrum in Figure 1. As before, we kept the ToA 
uncertainties and the sampling the same as for the actual 
observations of PSR J1713+0747 (note, here our data sets 
were not enlarged by 10,000 days as before). Our procedure 
is as follows: 

• Make 50 realizations of the noise (as in Section 3). 

• Measure the rms of the interpolation error between the 
different realisations for each sample. These are then aver- 
aged to give a single value. 

• This entire procedure is repeated 10 times. The mean 
and the rms of the mean are recorded. 

• The parameters of the model are changed (modifica- 
tions are made to a, f c or a\, in turn keeping the other 
two parameters fixed at their correct value) and the entire 
procedure repeated. 

The mean rms values are shown as a function of the 
error in the three model parameters as the three lines in 
Figure [6] The error bars (which are generally smaller than 
the symbol size) represent the rms of the mean. 

The largest discrepancies are seen with an incorrect de- 
termination of the corner frequency / c . The algorithm is very 
insensitive to errors in the spectral exponent, a. In all cases 
the discrepancies are significantly smaller than the mean 
rms of the white noise (0.69ps) and the mean rms of the 
red noise (0.42/is). We note that the parameter value that 
leads to the lowest mean rms in Figure [6] does not corre- 
spond exactly to the known model parameter. This is prob- 
ably because the red noise process is not stationary after the 
quadratic polynomial removal. We have also studied the im- 
pact of incorrect model parameters on extrapolating the tim- 
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ing residuals. For incorrect model parameters the duration 
for which the timing residuals can be extrapolated will be 
reduced when compared with the correct model parameters. 
The maximum extrapolation duration will depend upon the 
size of the error in the parameters, the white noise level and 
the requirements of the specific application. However, even 
with significant errors in the model parameters, the varia- 
tion in the discrepancy between the predicted residuals and 
the known simulated red noise will be smaller than varia- 
tions in different realisations of the red noise. The technique 
is therefore insensitive to errors in the model parameters and 
this will not affect our method for actual data sets. 



variance than the MLE, provided that the covariance ma- 
trix is correct. We have shown that in the pulsar timing 
context the MLE is not unduly sensitive to errors in the 
covariance matrix, so its performance must be close to op- 
timal. As data sets become longer our ability to model the 
timing noise will become even more important. Future ap- 
plications are already being planned that will require that 
the timing noise be extrapolated far into the future. It will 
be important to estimate the error incurred by such extrap- 
olations. The MLE estimator we have described provides a 
mechanism for doing so. 



5.3 Quadratic removal of timing noise 

The observed pulsar timing residuals depend on the method 
used to fit for the pulsar's pulse frequency and its first 
derivative. In F igure [7] (which is based on Figure 1 of 
IColes et al.ll201lh the upper section of the left-hand panel 
shows the timing residuals for PSR J1939+2134 obtained us- 
ing the Parkes radio telescope. The residuals are significantly 
red and, after fitting and removing a quadratic polynomial, 
take the form of a cubic polynomial. However, if this fit is 
carried out using the gen eralised least-squares-fitting pro- 
cedure l|Coles et al.ll201lT ) then the resulting residuals still 
exhibit a linear and quadratic term (this is shown in the 
right-panel of Figure [7}. We have used our algorithm to in- 
terpolate and extrapolate these residuals. The results are 
shown as the dotted line in the two panels of Figure [7] The 
difference between the model and the residuals are shown 
in the lower sections of each panel. In both cases the in- 
terpolation works well and converges on the predictions of 
the timing model. However, because of the different post-fit 
parameters, the timing models are different in the two cases. 

If this pulsar is being used for navigational purposes 
then measurements of absolute pulse arrival times are nec- 
essary. From the two panels of Figure [7] it may be expected 
that the two solutions and extrapolations will lead to dis- 
crepant measurements of the absolute pulse arrival times. 
In order to test this, we have determined pulse arrival times 
using both models shown in Figure [7J The difference be- 
tween the absolute time determinations using the two mod- 
els are shown in Figure [5] During the actual observations 
both models produce almost identical absolute pulse arrival 
times (with a discrepancy of less than 1 ns; shown in the 
lower panel). The extrapolated functions converge to the 
timing models; however, as these models are different, the 
predictions of the absolute arrival times diverge. The exact 
method in this case for measuring the pulsar's pulse fre- 
quency and its derivative is therefore not important when 
interpolating timing residuals to obtain absolute pulse ar- 
rival times, but does affect extrapolated arrival times. 



6 CONCLUSIONS 

In the presence of red timing noise we recommend that our 
new method be used to model, interpolate or extrapolate 
pulsar timing noise. This method is applicable even under 
extreme conditions such as very steep red spectra, very large 
gaps and/or highly variable ToA uncertainties. No other 
linear unbiased interpolator or extrapolator can have lower 
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Figure 7. Interpolation and extrapolation of the PSR J1939+2134 timing residuals (do tted line). The left -hand top panel shows the 
residuals after fitting for the pulse frequen cy and its first tim e derivative without using the lColes et al.l (|201 lh algorithm. The right-hand 
top panel shows the same after using the Coles et al. algorithm. The lower parts of each panel show the difference between the 

residuals and the interpolated function during the data span. 
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Figure 8. Comparison of two models for PSR J1939+2134. The 
top panels shows the difference between determinations of abso- 
lute arrival times given two models for PSR J1939+2134, shown 
in Figure 7. The bottom panel is a zoom-in of the upper panel 
during the region where actual data exist. 
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APPENDIX A: USAGE INSTRUCTIONS 

This algorithm has been implemented as part of the inter- 
polate plugin for tempo2. Using the plugin is straightfor- 
ward, but requires a model of the red timing noise. The noise 
is parameterised by / c , a and A according to: 



P(f) = 




The spectralModel plugin may be used to model the red 
noise component of the timing residuals and determine these 
parameters. 

Once these parameters are determined the interpo- 
late plugin can be used as follows: 

> tempo2 -gr interpolate -f mypar.par mytim.tim 

-a A -fc fc -alpha alpha 

By default the software produces a new parameter file that 
includes a set of parameters that give the interpolated func- 
tion at specified times. These parameters are defined using 
the IFUNC keyword. It is also possible to provide differ- 
ent sampling for purposes of e.g., extrapolation using the 
-xl, -x2 and -dx command-line options. These correspond 
to the start time, the last time and the step size for the 
extrapolation, respectively. 

Figures 4 and 8 were produced using the plotInterp 
which allows the user to visualise the interpolation. Typical 
usage is: 

> tempo2 -gr plotInterp -f mypar.par mytim.tim 

-xextend 5000 
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The xextend command line option extrapolates the predic- 
tion for 5000 d before the first observation and 5000 d after 
the last observation. 
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